Index

Library setup

suppressMessages(
  suppressWarnings(
    c(library(data.table),
      library(gdata),
      library(ggplot2),
      library(RColorBrewer),
      library(readr),
      library(tidyverse),
      library(plotly),
      library(directlabels),
      library(gplots),
      library(knitr)
      )))
##   [1] "data.table"   "stats"        "graphics"     "grDevices"    "utils"       
##   [6] "datasets"     "methods"      "base"         "gdata"        "data.table"  
##  [11] "stats"        "graphics"     "grDevices"    "utils"        "datasets"    
##  [16] "methods"      "base"         "ggplot2"      "gdata"        "data.table"  
##  [21] "stats"        "graphics"     "grDevices"    "utils"        "datasets"    
##  [26] "methods"      "base"         "RColorBrewer" "ggplot2"      "gdata"       
##  [31] "data.table"   "stats"        "graphics"     "grDevices"    "utils"       
##  [36] "datasets"     "methods"      "base"         "readr"        "RColorBrewer"
##  [41] "ggplot2"      "gdata"        "data.table"   "stats"        "graphics"    
##  [46] "grDevices"    "utils"        "datasets"     "methods"      "base"        
##  [51] "forcats"      "stringr"      "dplyr"        "purrr"        "tidyr"       
##  [56] "tibble"       "tidyverse"    "readr"        "RColorBrewer" "ggplot2"     
##  [61] "gdata"        "data.table"   "stats"        "graphics"     "grDevices"   
##  [66] "utils"        "datasets"     "methods"      "base"         "plotly"      
##  [71] "forcats"      "stringr"      "dplyr"        "purrr"        "tidyr"       
##  [76] "tibble"       "tidyverse"    "readr"        "RColorBrewer" "ggplot2"     
##  [81] "gdata"        "data.table"   "stats"        "graphics"     "grDevices"   
##  [86] "utils"        "datasets"     "methods"      "base"         "directlabels"
##  [91] "plotly"       "forcats"      "stringr"      "dplyr"        "purrr"       
##  [96] "tidyr"        "tibble"       "tidyverse"    "readr"        "RColorBrewer"
## [101] "ggplot2"      "gdata"        "data.table"   "stats"        "graphics"    
## [106] "grDevices"    "utils"        "datasets"     "methods"      "base"        
## [111] "gplots"       "directlabels" "plotly"       "forcats"      "stringr"     
## [116] "dplyr"        "purrr"        "tidyr"        "tibble"       "tidyverse"   
## [121] "readr"        "RColorBrewer" "ggplot2"      "gdata"        "data.table"  
## [126] "stats"        "graphics"     "grDevices"    "utils"        "datasets"    
## [131] "methods"      "base"         "knitr"        "gplots"       "directlabels"
## [136] "plotly"       "forcats"      "stringr"      "dplyr"        "purrr"       
## [141] "tidyr"        "tibble"       "tidyverse"    "readr"        "RColorBrewer"
## [146] "ggplot2"      "gdata"        "data.table"   "stats"        "graphics"    
## [151] "grDevices"    "utils"        "datasets"     "methods"      "base"

ID conversion from human gene symbol to mouse Ensembl ID

This conversion is done using ID database from MGI to be consistent with Doug Brubaker’s analysis using human RNA-Seq dataset.

Import the dataset containing Gene ID information from MGI

homolog_id <- read_csv("Human_2_mouse_symbol.csv")
## Parsed with column specification:
## cols(
##   `Human Marker Symbol` = col_character(),
##   `Human Entrez Gene ID` = col_double(),
##   `Mouse Marker Symbol` = col_character(),
##   `MGI Marker Accession ID` = col_character()
## )
ensembl_id <- read_csv("mouse_symbol_Ensembl.csv")
## Parsed with column specification:
## cols(
##   `MGI Accession ID` = col_character(),
##   `Marker Symbol` = col_character(),
##   `EntrezGene ID` = col_character(),
##   `Ensembl Gene ID` = col_character(),
##   `HomoloGene ID` = col_character()
## )
# load the PC dataset
cCD_PC <- read_csv("cCD PC.csv")
## Parsed with column specification:
## cols(
##   Row = col_character(),
##   cCD_PC1 = col_double(),
##   cCD_PC2 = col_double(),
##   cCD_PC3 = col_double(),
##   cCD_PC4 = col_double(),
##   cCD_PC5 = col_double(),
##   cCD_PC6 = col_double(),
##   cCD_PC7 = col_double(),
##   cCD_PC8 = col_double(),
##   cCD_PC9 = col_double(),
##   cCD_PC10 = col_double(),
##   cCD_PC11 = col_double(),
##   cCD_PC12 = col_double(),
##   cCD_PC13 = col_double(),
##   cCD_PC14 = col_double(),
##   cCD_PC15 = col_double(),
##   cCD_PC16 = col_double()
## )
iCD_PC <- read_csv("iCD PC.csv")
## Parsed with column specification:
## cols(
##   Row = col_character(),
##   iCD_PC1 = col_double(),
##   iCD_PC2 = col_double(),
##   iCD_PC3 = col_double(),
##   iCD_PC4 = col_double(),
##   iCD_PC5 = col_double(),
##   iCD_PC6 = col_double(),
##   iCD_PC7 = col_double(),
##   iCD_PC8 = col_double(),
##   iCD_PC9 = col_double(),
##   iCD_PC10 = col_double(),
##   iCD_PC11 = col_double(),
##   iCD_PC12 = col_double(),
##   iCD_PC13 = col_double(),
##   iCD_PC14 = col_double(),
##   iCD_PC15 = col_double(),
##   iCD_PC16 = col_double()
## )
human_id <- as.data.frame(cCD_PC$Row)
human_id <- as_tibble(human_id)
colnames(human_id) <- "human_gene_symbol"

id_conversion <- inner_join(human_id, homolog_id, by = c("human_gene_symbol" = "Human Marker Symbol"))
## Warning: Column `human_gene_symbol`/`Human Marker Symbol` joining factor and
## character vector, coercing into character vector
non_duplicates_idx <- which(duplicated(id_conversion$human_gene_symbol) == FALSE)
id_conversion <- id_conversion[non_duplicates_idx, ]

id_conversion <- inner_join(id_conversion, ensembl_id, by = c("Mouse Marker Symbol" = "Marker Symbol"))
non_duplicates_idx <- which(duplicated(id_conversion$`Mouse Marker Symbol`) == FALSE)
id_conversion <- id_conversion[non_duplicates_idx, ]

id_conversion <- id_conversion[,c(1,7)]

cCD_PC <- inner_join(cCD_PC, id_conversion, by = c("Row" = "human_gene_symbol"))
iCD_PC <- inner_join(iCD_PC, id_conversion, by = c("Row" = "human_gene_symbol"))

Heatmap plotting based on cCD PC

For cCD PCs, PC1 is selected because it seperates diseases statuses. PC10 and PC14 are selected because they seperate TNF from TCT mice.

The goal of this is to examine how many of the genes that contribute to the seperation of TNF and TCT model based on human PCs are targeted by MK2 inhibitor.

For each PC, the top 325 contributing genes are selected for plotting the heatmap (10% of the list). Top contributing genes are genes with the largest absolute loading coefficients.

n = 10
threshold_p1 <- quantile(abs(cCD_PC$cCD_PC1), prob=1-n/100)
threshold_p10 <- quantile(abs(cCD_PC$cCD_PC10), prob=1-n/100)
threshold_p14 <- quantile(abs(cCD_PC$cCD_PC14), prob=1-n/100)

ensembl_p1 <- cCD_PC$`Ensembl Gene ID`[abs(cCD_PC$cCD_PC1) > threshold_p1]
ensembl_p10 <- cCD_PC$`Ensembl Gene ID`[abs(cCD_PC$cCD_PC10) > threshold_p10]
ensembl_p14 <- cCD_PC$`Ensembl Gene ID`[abs(cCD_PC$cCD_PC14) > threshold_p14]

cCD_PC1

I will use I_V vs NI_V data to plot this heatmap since this is the PC that separates diseases status.

TCT_iv_niv_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/VST_I_V vs NI_V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   NI_V_1 = col_double(),
##   NI_V_2 = col_double(),
##   NI_V_4 = col_double(),
##   NI_V_6 = col_double(),
##   NI_V_7 = col_double(),
##   I_V_1 = col_double(),
##   I_V_2 = col_double(),
##   I_V_3 = col_double(),
##   I_V_4 = col_double(),
##   I_V_5 = col_double(),
##   I_V_6 = col_double()
## )
TNF_iv_niv_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/VST_I_V vs NI_V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   NI_V_1 = col_double(),
##   NI_V_2 = col_double(),
##   NI_V_3 = col_double(),
##   NI_V_4 = col_double(),
##   NI_V_5 = col_double(),
##   I_V_1 = col_double(),
##   I_V_2 = col_double(),
##   I_V_3 = col_double(),
##   I_V_4 = col_double(),
##   I_V_5 = col_double()
## )
TCT_iv_niv_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Differential Analysis_I_V vs NI_V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )
TNF_iv_niv_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Differential Analysis_I_V vs NI_V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )

TCT model

All genes in the list of 325 from PC1

suppressMessages(library(mosaic))

pc1_genes <- TCT_iv_niv_vst[TCT_iv_niv_vst$X1 %in% ensembl_p1,]
pc1_genes <- as.data.frame(pc1_genes)

for (i in 1:dim(pc1_genes)[1]) {
  pc1_genes[i,2:12] <- zscore(as.numeric(pc1_genes[i,2:12]))
}

my_palette <- colorRampPalette(c("blue", "white", "red"))(128)
heatmap_matrix <- as.matrix(pc1_genes[,2:12])

png('TCT_I_V vs NI_V in cCD PC1 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)

heatmap.2(heatmap_matrix,
          main = "TCT I_V vs NI_V\nin cCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_I_V vs NI_V in cCD PC1 Genes.pdf',
    width = 6,
    height = 10)

heatmap.2(heatmap_matrix,
          main = "TCT I_V vs NI_V\nin cCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics("TCT_I_V vs NI_V in cCD PC1 Genes.png")

DEG from PC1

pc1_dif <- TCT_iv_niv_dif[TCT_iv_niv_dif$X1 %in% ensembl_p1,]
pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.05  & abs(pc1_dif$log2FoldChange) > 0.25)
pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))


heatmap_matrix <- as.matrix(pc1_genes_DEG[,2:12])

png('TCT_I_V vs NI_V DEG in cCD PC1 Genes.png',
    width = 800,
    height = 600,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix,
          main = "TCT_I_V vs NI_V DEG\nin cCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,3),
          col=my_palette,
          cexCol = 1,
          margins = c(4,6),
          trace = "none",
          dendrogram = "row",
          labRow = pc1_genes_DEG$`Marker Symbol`,
          keysize = 2,
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_I_V vs NI_V DEG in cCD PC1 Genes.pdf',
    width = 8,
    height = 6)
heatmap.2(heatmap_matrix,
          main = "TCT_I_V vs NI_V DEG\nin cCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,3),
          col=my_palette,
          cexCol = 1,
          margins = c(4,6),
          trace = "none",
          dendrogram = "row",
          labRow = pc1_genes_DEG$`Marker Symbol`,
          keysize = 2,
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics("TCT_I_V vs NI_V DEG in cCD PC1 Genes.png")

6/324 genes are deemed significant.

TNF model

All genes in the list of 325 from PC1

pc1_genes <- TNF_iv_niv_vst[TNF_iv_niv_vst$X1 %in% ensembl_p1,]
pc1_genes <- as.data.frame(pc1_genes)


for (i in 1:dim(pc1_genes)[1]) {
  pc1_genes[i,2:11] <- zscore(as.numeric(pc1_genes[i,2:11]))
}

heatmap_matrix <- as.matrix(pc1_genes[,2:11])

png('TNF_I_V vs NI_V in cCD PC1 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix,
          main = "TNF I_V vs NI_V\nin cCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_I_V vs NI_V in cCD PC1 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix,
          main = "TNF I_V vs NI_V\nin cCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics("TNF_I_V vs NI_V in cCD PC1 Genes.png")

DEG from PC1

pc1_dif <- TNF_iv_niv_dif[TNF_iv_niv_dif$X1 %in% ensembl_p1,]
pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.05 & abs(pc1_dif$log2FoldChange) > 0.25)
pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))


heatmap_matrix <- as.matrix(pc1_genes_DEG[,2:11])

png('TNF_I_V vs NI_V DEG in cCD PC1 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix,
          main = "TNF\nI_V vs NI_V DEG\ncCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(4,6),
          trace = "none",
          dendrogram = "row",
          labRow = pc1_genes_DEG$`Marker Symbol`,
          keysize = 2,
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_I_V vs NI_V DEG in cCD PC1 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix,
          main = "TNF\nI_V vs NI_V DEG\ncCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(4,6),
          trace = "none",
          dendrogram = "row",
          labRow = pc1_genes_DEG$`Marker Symbol`,
          keysize = 2,
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_I_V vs NI_V DEG in cCD PC1 Genes.png')

74/324 genes are deemed significant.

cCD_PC10 and cCD_PC14

I_MKI vs I_V datasets will be used here.

TCT_imki_iv_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/VST_I_MKI vs I_V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   I_V_1 = col_double(),
##   I_V_2 = col_double(),
##   I_V_3 = col_double(),
##   I_V_4 = col_double(),
##   I_V_5 = col_double(),
##   I_V_6 = col_double(),
##   I_MKI_1 = col_double(),
##   I_MKI_2 = col_double(),
##   I_MKI_3 = col_double(),
##   I_MKI_4 = col_double(),
##   I_MKI_5 = col_double(),
##   I_MKI_6 = col_double()
## )
TCT_nimki_niv_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/VST_NI_MKI vs NI_V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   NI_V_1 = col_double(),
##   NI_V_2 = col_double(),
##   NI_V_4 = col_double(),
##   NI_V_6 = col_double(),
##   NI_V_7 = col_double(),
##   NI_MKI_1 = col_double(),
##   NI_MKI_2 = col_double(),
##   NI_MKI_3 = col_double(),
##   NI_MKI_4 = col_double(),
##   NI_MKI_5 = col_double(),
##   NI_MKI_6 = col_double()
## )
TCT_combined_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/VST_MKI vs V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   X1 = col_character()
## )
## See spec(...) for full column specifications.
TNF_imki_iv_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/VST_I_MKI vs I_V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   I_V_1 = col_double(),
##   I_V_2 = col_double(),
##   I_V_3 = col_double(),
##   I_V_4 = col_double(),
##   I_V_5 = col_double(),
##   I_MKI_1 = col_double(),
##   I_MKI_2 = col_double(),
##   I_MKI_3 = col_double(),
##   I_MKI_4 = col_double(),
##   I_MKI_5 = col_double(),
##   I_MKI_6 = col_double()
## )
TNF_nimki_niv_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/VST_NI_MKI vs NI_V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   NI_V_1 = col_double(),
##   NI_V_2 = col_double(),
##   NI_V_3 = col_double(),
##   NI_V_4 = col_double(),
##   NI_V_5 = col_double(),
##   NI_MKI_1 = col_double(),
##   NI_MKI_2 = col_double(),
##   NI_MKI_3 = col_double(),
##   NI_MKI_4 = col_double(),
##   NI_MKI_5 = col_double()
## )
TNF_combined_vst <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/VST_MKI vs V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   X1 = col_character()
## )
## See spec(...) for full column specifications.
TCT_imki_iv_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Differential Analysis_I_MKI vs I_V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )
TCT_nimki_niv_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Differential Analysis_NI_MKI vs NI_V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )
TCT_combined_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TCT/Differential Analysis_MKI vs V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )
TNF_imki_iv_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Differential Analysis_I_MKI vs I_V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )
TNF_nimki_niv_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Differential Analysis_NI_MKI vs NI_V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )
TNF_combined_dif <- read_csv("~/OneDrive - Harvard University/Haigis Lab/Projects/Colloboration Projects/IBD RNA-Seq analysis/Analysis/TNF/Differential Analysis_MKI vs V.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_character(),
##   baseMean = col_double(),
##   log2FoldChange = col_double(),
##   lfcSE = col_double(),
##   stat = col_double(),
##   pvalue = col_double(),
##   padj = col_double()
## )

TCT model stratifying based on inflammation status

All genes in the list of 325 from PC10 and PC14

pc10_genes <- as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p10,])
pc14_genes <- as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p14,])


for (i in 1:dim(pc10_genes)[1]) {
  pc10_genes[i,2:13] <- zscore(as.numeric(pc10_genes[i,2:13]))
}

for (i in 1:dim(pc14_genes)[1]) {
  pc14_genes[i,2:13] <- zscore(as.numeric(pc14_genes[i,2:13]))
}

heatmap_matrix_pc10 <- as.matrix(pc10_genes[,2:13])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:13])


png('TCT_I_MKI vs I_V in cCD PC10 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc10,
          main = "TCT I_MKI vs I_V\nin cCD PC10 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_I_MKI vs I_V in cCD PC10 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc10,
          main = "TCT I_MKI vs I_V\nin cCD PC10 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TCT_I_MKI vs I_V in cCD PC10 Genes.png')

png('TCT_I_MKI vs I_V in cCD PC14 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
          main = "TCT I_MKI vs I_V\nin cCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_I_MKI vs I_V in cCD PC14 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc14,
          main = "TCT I_MKI vs I_V\nin cCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TCT_I_MKI vs I_V in cCD PC14 Genes.png')

DEG from PC10 and PC14. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.

pc10_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p10,]
pc14_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p14,]

pc10_DEG <- subset(pc10_dif, pc10_dif$padj < 0.1  & abs(pc10_dif$log2FoldChange > 0.1))
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1  & abs(pc14_dif$log2FoldChange > 0.1))

pc10_genes_DEG <- pc10_genes[pc10_genes$X1 %in% pc10_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]

pc10_genes_DEG <- inner_join(pc10_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))

print(paste("The number of MK2 inhibitor targets from TCT model from cCD PC10:", dim(pc10_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TCT model from cCD PC10: 0"
print(paste("The number of MK2 inhibitor targets from TCT model from cCD PC14:", dim(pc14_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TCT model from cCD PC14: 1"

We cannot make heatmap with only 1 item as input. So I will just print out the item here.

pc14_DEG
## # A tibble: 1 x 7
##   X1                 baseMean log2FoldChange lfcSE  stat     pvalue     padj
##   <chr>                 <dbl>          <dbl> <dbl> <dbl>      <dbl>    <dbl>
## 1 ENSMUSG00000050022     55.3           1.34 0.263  4.82 0.00000141 0.000974

I would also like to plot the Loading plot for this gene in PC10 and PC14 space.

loading_cord <- cCD_PC[cCD_PC$`Ensembl Gene ID` %in% pc14_genes_DEG$X1,]
loading_cord <- loading_cord[, c('cCD_PC10', 'cCD_PC14', 'Ensembl Gene ID')]
loading_cord <- inner_join(loading_cord, ensembl_id[,c(2,4)], by = c("Ensembl Gene ID"="Ensembl Gene ID"))
loading_cord$`Ensembl Gene ID` = NULL
loading_cord[2,] <- loading_cord[1,]
loading_cord <- as.data.frame(loading_cord)
loading_cord[1,] <- c(0,0,NA)
loading_cord$cCD_PC10 <- as.numeric(loading_cord$cCD_PC10)
loading_cord$cCD_PC14 <- as.numeric(loading_cord$cCD_PC14)


ggplot(loading_cord, aes(x=cCD_PC14, y=cCD_PC10)) + 
  geom_line(color = "red",
            arrow = arrow()) +
  xlim(-0.1, 0.1) +
  ylim(-0.1, 0.1) +
  ggtitle("Loading plot for TCT MK2i target in cCD PC10 and PC14 space") +
  geom_dl(aes(label = loading_cord$`Marker Symbol`), method = list(dl.trans(x = x + 0.2), "last.points", cex = 0.8))
## Warning: Use of `loading_cord$`Marker Symbol`` is discouraged. Use `Marker
## Symbol` instead.

TCT model disregarding inflammation status

All genes in the list of 325 from PC10 and PC14

pc10_genes <- as.data.frame(TCT_combined_vst[TCT_combined_vst$X1 %in% ensembl_p10,])
pc14_genes <- as.data.frame(TCT_combined_vst[TCT_combined_vst$X1 %in% ensembl_p14,])


for (i in 1:dim(pc10_genes)[1]) {
  pc10_genes[i,2:26] <- zscore(as.numeric(pc10_genes[i,2:26]))
}

for (i in 1:dim(pc14_genes)[1]) {
  pc14_genes[i,2:26] <- zscore(as.numeric(pc14_genes[i,2:26]))
}

heatmap_matrix_pc10 <- as.matrix(pc10_genes[,2:26])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:26])


png('TCT_MKI vs V in cCD PC10 Genes.png',
    width = 800,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc10,
          main = "TCT MKI vs V\nin cCD PC10 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_MKI vs V in cCD PC10 Genes.pdf',
    width = 8,
    height = 10)
heatmap.2(heatmap_matrix_pc10,
          main = "TCT MKI vs V\nin cCD PC10 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TCT_MKI vs V in cCD PC10 Genes.png')

png('TCT_MKI vs V in cCD PC14 Genes.png',
    width = 800,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
          main = "TCT MKI vs V\nin cCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_MKI vs V in cCD PC14 Genes.pdf',
    width = 8,
    height = 10)
heatmap.2(heatmap_matrix_pc14,
          main = "TCT MKI vs V\nin cCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TCT_MKI vs V in cCD PC14 Genes.png')

DEG from PC10 and PC14. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.

pc10_dif <- TCT_combined_dif[TCT_combined_dif$X1 %in% ensembl_p10,]
pc14_dif <- TCT_combined_dif[TCT_combined_dif$X1 %in% ensembl_p14,]

pc10_DEG <- subset(pc10_dif, pc10_dif$padj < 0.1 & abs(pc10_dif$log2FoldChange) > 0.1)
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1 & abs(pc14_dif$log2FoldChange) > 0.1)

pc10_genes_DEG <- pc10_genes[pc10_genes$X1 %in% pc10_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]

pc10_genes_DEG <- inner_join(pc10_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))

print(paste("The number of MK2 inhibitor targets from combined TCT model from cCD PC10:", dim(pc10_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from combined TCT model from cCD PC10: 0"
print(paste("The number of MK2 inhibitor targets from combined TCT model from cCD PC14:", dim(pc14_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from combined TCT model from cCD PC14: 0"

TNF model stratifying based on inflammation status

All genes in the list of 325 from PC10 and PC14

pc10_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p10,])
pc14_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p14,])


for (i in 1:dim(pc10_genes)[1]) {
  pc10_genes[i,2:12] <- zscore(as.numeric(pc10_genes[i,2:12]))
}

for (i in 1:dim(pc14_genes)[1]) {
  pc14_genes[i,2:12] <- zscore(as.numeric(pc14_genes[i,2:12]))
}

heatmap_matrix_pc10 <- as.matrix(pc10_genes[,2:12])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:12])


png('TNF_I_MKI vs I_V in cCD PC10 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc10,
          main = "TNF I_MKI vs I_V\nin cCD PC10 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_I_MKI vs I_V in cCD PC10 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc10,
          main = "TNF I_MKI vs I_V\nin cCD PC10 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_I_MKI vs I_V in cCD PC10 Genes.png')

png('TNF_I_MKI vs I_V in cCD PC14 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
          main = "TNF I_MKI vs I_V\nin cCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_I_MKI vs I_V in cCD PC14 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc14,
          main = "TNF I_MKI vs I_V\nin cCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_I_MKI vs I_V in cCD PC14 Genes.png')

DEG from PC10 and PC14. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.

pc10_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p10,]
pc14_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p14,]

pc10_DEG <- subset(pc10_dif, pc10_dif$padj < 0.1  & abs(pc10_dif$log2FoldChange > 0.1))
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1  & abs(pc14_dif$log2FoldChange > 0.1))

pc10_genes_DEG <- pc10_genes[pc10_genes$X1 %in% pc10_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]

pc10_genes_DEG <- inner_join(pc10_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))

print(paste("The number of MK2 inhibitor targets from TNF model from cCD PC10:", dim(pc10_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from cCD PC10: 0"
print(paste("The number of MK2 inhibitor targets from TNF model from cCD PC14:", dim(pc14_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from cCD PC14: 0"

TNF model disregarding inflammation status

All genes in the list of 325 from PC10 and PC14

pc10_genes <- as.data.frame(TNF_combined_vst[TNF_combined_vst$X1 %in% ensembl_p10,])
pc14_genes <- as.data.frame(TNF_combined_vst[TNF_combined_vst$X1 %in% ensembl_p14,])


for (i in 1:dim(pc10_genes)[1]) {
  pc10_genes[i,2:22] <- zscore(as.numeric(pc10_genes[i,2:22]))
}

for (i in 1:dim(pc14_genes)[1]) {
  pc14_genes[i,2:22] <- zscore(as.numeric(pc14_genes[i,2:22]))
}

heatmap_matrix_pc10 <- as.matrix(pc10_genes[,2:22])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:22])


png('TNF_MKI vs V in cCD PC10 Genes.png',
    width = 800,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc10,
          main = "TNF MKI vs V\nin cCD PC10 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_MKI vs V in cCD PC10 Genes.pdf',
    width = 8,
    height = 10)
heatmap.2(heatmap_matrix_pc10,
          main = "TNF MKI vs V\nin cCD PC10 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_MKI vs V in cCD PC10 Genes.png')

png('TNF_MKI vs V in cCD PC14 Genes.png',
    width = 800,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
          main = "TNF MKI vs V\nin cCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_MKI vs V in cCD PC14 Genes.pdf',
    width = 8,
    height = 10)
heatmap.2(heatmap_matrix_pc14,
          main = "TNF MKI vs V\nin cCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_MKI vs V in cCD PC14 Genes.png')

DEG from PC10 and PC14. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.

pc10_dif <- TNF_combined_dif[TNF_combined_dif$X1 %in% ensembl_p10,]
pc14_dif <- TNF_combined_dif[TNF_combined_dif$X1 %in% ensembl_p14,]

pc10_DEG <- subset(pc10_dif, pc10_dif$padj < 0.1 & abs(pc10_dif$log2FoldChange) > 0.1)
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1 & abs(pc14_dif$log2FoldChange) > 0.1)

pc10_genes_DEG <- pc10_genes[pc10_genes$X1 %in% pc10_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]

pc10_genes_DEG <- inner_join(pc10_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))

print(paste("The number of MK2 inhibitor targets from TNF model from cCD PC10:", dim(pc10_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from cCD PC10: 1"
print(paste("The number of MK2 inhibitor targets from TNF model from cCD PC14:", dim(pc14_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from cCD PC14: 0"

We cannot make heatmap with only 1 item as input. So I will just print out the item here.

pc10_DEG
## # A tibble: 1 x 7
##   X1                 baseMean log2FoldChange lfcSE  stat      pvalue     padj
##   <chr>                 <dbl>          <dbl> <dbl> <dbl>       <dbl>    <dbl>
## 1 ENSMUSG00000028358     102.          0.582 0.175  4.94 0.000000781 0.000874

I would also like to plot the Loading plot for this gene in PC10 and PC14 space.

loading_cord <- cCD_PC[cCD_PC$`Ensembl Gene ID` %in% pc10_genes_DEG$X1,]
loading_cord <- loading_cord[, c('cCD_PC10', 'cCD_PC14', 'Ensembl Gene ID')]
loading_cord <- inner_join(loading_cord, ensembl_id[,c(2,4)], by = c("Ensembl Gene ID"="Ensembl Gene ID"))
loading_cord$`Ensembl Gene ID` = NULL
loading_cord[2,] <- loading_cord[1,]
loading_cord <- as.data.frame(loading_cord)
loading_cord[1,] <- c(0,0,NA)
loading_cord$cCD_PC10 <- as.numeric(loading_cord$cCD_PC10)
loading_cord$cCD_PC14 <- as.numeric(loading_cord$cCD_PC14)


ggplot(loading_cord, aes(x=cCD_PC14, y=cCD_PC10)) + 
  geom_line(color = "red",
            arrow = arrow(ends = "first")) +
  xlim(-0.1, 0.1) +
  ylim(-0.1, 0.1) +
  ggtitle("Loading plot for combined TNF MK2i target in cCD PC10 and PC14 space") +
  geom_dl(aes(label = loading_cord$`Marker Symbol`), method = list(dl.trans(x = x + 0.2), "last.points", cex = 0.8))
## Warning: Use of `loading_cord$`Marker Symbol`` is discouraged. Use `Marker
## Symbol` instead.

Heatmap plotting based on iCD PC

For iCD PCs, PC1, PC14 and PC16 are selected because they seperate TNF from TCT mice.

The goal of this is to examine how many of the genes that contribute to the seperation of TNF and TCT model based on human PCs are targeted by MK2 inhibitor.

For each PC, the top 325 contributing genes are selected for plotting the heatmap (10% of the list). Top contributing genes are genes with the largest absolute loading coefficients.

threshold_p1 <- quantile(abs(iCD_PC$iCD_PC1), prob=1-n/100)
threshold_p14 <- quantile(abs(iCD_PC$iCD_PC14), prob=1-n/100)
threshold_p16 <- quantile(abs(iCD_PC$iCD_PC16), prob=1-n/100)

ensembl_p1 <- iCD_PC$`Ensembl Gene ID`[abs(iCD_PC$iCD_PC1) > threshold_p1]
ensembl_p14 <- iCD_PC$`Ensembl Gene ID`[abs(iCD_PC$iCD_PC14) > threshold_p14]
ensembl_p16 <- iCD_PC$`Ensembl Gene ID`[abs(iCD_PC$iCD_PC16) > threshold_p16]

I_MKI vs I_V datasets will be used here for the analysis.

TCT model stratifying based on inflammation status

All genes in the list of 325 from PC1, PC14 and PC16

pc1_genes <- as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p1,])
pc14_genes <- as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p14,])
pc16_genes <-as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p16,])


for (i in 1:dim(pc1_genes)[1]) {
  pc1_genes[i,2:13] <- zscore(as.numeric(pc1_genes[i,2:13]))
}

for (i in 1:dim(pc14_genes)[1]) {
  pc14_genes[i,2:13] <- zscore(as.numeric(pc14_genes[i,2:13]))
}

for (i in 1:dim(pc16_genes)[1]) {
  pc16_genes[i,2:13] <- zscore(as.numeric(pc16_genes[i,2:13]))
}


heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:13])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:13])
heatmap_matrix_pc16 <- as.matrix(pc16_genes[,2:13])


png('TCT_I_MKI vs I_V in iCD PC1 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
          main = "TCT I_MKI vs I_V\nin iCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_I_MKI vs I_V in iCD PC1 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc1,
          main = "TCT I_MKI vs I_V\nin iCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TCT_I_MKI vs I_V in iCD PC1 Genes.png')

png('TCT_I_MKI vs I_V in iCD PC14 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
          main = "TCT I_MKI vs I_V\nin iCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_I_MKI vs I_V in iCD PC14 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc14,
          main = "TCT I_MKI vs I_V\nin iCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TCT_I_MKI vs I_V in iCD PC14 Genes.png')

png('TCT_I_MKI vs I_V in iCD PC16 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc16,
          main = "TCT I_MKI vs I_V\nin iCD PC16 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_I_MKI vs I_V in iCD PC16 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc16,
          main = "TCT I_MKI vs I_V\nin iCD PC16 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TCT_I_MKI vs I_V in iCD PC16 Genes.png')

DEG from PC1, PC14 and PC16. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.

pc1_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p1,]
pc14_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p14,]
pc16_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p16,]


pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1  & abs(pc1_dif$log2FoldChange > 0.1))
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1  & abs(pc14_dif$log2FoldChange > 0.1))
pc16_DEG <- subset(pc16_dif, pc16_dif$padj < 0.1  & abs(pc16_dif$log2FoldChange > 0.1))


pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]
pc16_genes_DEG <- pc16_genes[pc16_genes$X1 %in% pc16_DEG$X1,]


pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc16_genes_DEG <- inner_join(pc16_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))


print(paste("The number of MK2 inhibitor targets from TCT model from iCD PC1:", dim(pc1_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TCT model from iCD PC1: 0"
print(paste("The number of MK2 inhibitor targets from TCT model from iCD PC14:", dim(pc14_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TCT model from iCD PC14: 0"
print(paste("The number of MK2 inhibitor targets from TCT model from iCD PC16:", dim(pc16_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TCT model from iCD PC16: 0"

TCT model disregarding inflammation status

All genes in the list of 325 from PC1, PC14 and PC16

pc1_genes <- as.data.frame(TCT_combined_vst[TCT_combined_vst$X1 %in% ensembl_p1,])
pc14_genes <- as.data.frame(TCT_combined_vst[TCT_combined_vst$X1 %in% ensembl_p14,])
pc16_genes <-as.data.frame(TCT_combined_vst[TCT_combined_vst$X1 %in% ensembl_p16,])


for (i in 1:dim(pc1_genes)[1]) {
  pc1_genes[i,2:26] <- zscore(as.numeric(pc1_genes[i,2:26]))
}

for (i in 1:dim(pc14_genes)[1]) {
  pc14_genes[i,2:26] <- zscore(as.numeric(pc14_genes[i,2:26]))
}

for (i in 1:dim(pc16_genes)[1]) {
  pc16_genes[i,2:26] <- zscore(as.numeric(pc16_genes[i,2:26]))
}


heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:26])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:26])
heatmap_matrix_pc16 <- as.matrix(pc16_genes[,2:26])


png('TCT_MKI vs V in iCD PC1 Genes.png',
    width = 800,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
          main = "TCT MKI vs V\nin iCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_MKI vs V in iCD PC1 Genes.pdf',
    width = 8,
    height = 10)
heatmap.2(heatmap_matrix_pc1,
          main = "TCT MKI vs V\nin iCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics("TCT_MKI vs V in iCD PC1 Genes.png")

png('TCT_MKI vs V in iCD PC14 Genes.png',
    width = 800,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
          main = "TCT MKI vs V\nin iCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_MKI vs V in iCD PC14 Genes.pdf',
    width = 8,
    height = 10)
heatmap.2(heatmap_matrix_pc14,
          main = "TCT MKI vs V\nin iCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics("TCT_MKI vs V in iCD PC14 Genes.png")

png('TCT_MKI vs V in iCD PC16 Genes.png',
    width = 800,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc16,
          main = "TCT MKI vs V\nin iCD PC16 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_MKI vs V in iCD PC16 Genes.pdf',
    width = 8,
    height = 10)
heatmap.2(heatmap_matrix_pc16,
          main = "TCT MKI vs V\nin iCD PC16 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics("TCT_MKI vs V in iCD PC16 Genes.png")

DEG from PC1, PC14 and PC16. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.

pc1_dif <- TCT_combined_dif[TCT_combined_dif$X1 %in% ensembl_p1,]
pc14_dif <- TCT_combined_dif[TCT_combined_dif$X1 %in% ensembl_p14,]
pc16_dif <- TCT_combined_dif[TCT_combined_dif$X1 %in% ensembl_p16,]


pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1  & abs(pc1_dif$log2FoldChange > 0.1))
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1  & abs(pc14_dif$log2FoldChange > 0.1))
pc16_DEG <- subset(pc16_dif, pc16_dif$padj < 0.1  & abs(pc16_dif$log2FoldChange > 0.1))


pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]
pc16_genes_DEG <- pc16_genes[pc16_genes$X1 %in% pc16_DEG$X1,]


pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc16_genes_DEG <- inner_join(pc16_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))


print(paste("The number of MK2 inhibitor targets from TCT model from iCD PC1:", dim(pc1_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TCT model from iCD PC1: 0"
print(paste("The number of MK2 inhibitor targets from TCT model from iCD PC14:", dim(pc14_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TCT model from iCD PC14: 0"
print(paste("The number of MK2 inhibitor targets from TCT model from iCD PC16:", dim(pc16_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TCT model from iCD PC16: 0"

TNF model stratifying based on inflammation status

All genes in the list of 325 from PC1, PC14 and PC16

pc1_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p1,])
pc14_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p14,])
pc16_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p16,])


for (i in 1:dim(pc1_genes)[1]) {
  pc1_genes[i,2:12] <- zscore(as.numeric(pc1_genes[i,2:12]))
}

for (i in 1:dim(pc14_genes)[1]) {
  pc14_genes[i,2:12] <- zscore(as.numeric(pc14_genes[i,2:12]))
}

for (i in 1:dim(pc16_genes)[1]) {
  pc16_genes[i,2:12] <- zscore(as.numeric(pc16_genes[i,2:12]))
}


heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:12])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:12])
heatmap_matrix_pc16 <- as.matrix(pc16_genes[,2:12])



png('TNF_I_MKI vs I_V in iCD PC1 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
          main = "TNF I_MKI vs I_V\nin iCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_I_MKI vs I_V in iCD PC1 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc1,
          main = "TNF I_MKI vs I_V\nin iCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_I_MKI vs I_V in iCD PC1 Genes.png')

png('TNF_I_MKI vs I_V in iCD PC14 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
          main = "TNF I_MKI vs I_V\nin iCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_I_MKI vs I_V in iCD PC14 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc14,
          main = "TNF I_MKI vs I_V\nin iCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_I_MKI vs I_V in iCD PC14 Genes.png')

png('TNF_I_MKI vs I_V in iCD PC16 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc16,
          main = "TNF I_MKI vs I_V\nin iCD PC16 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_I_MKI vs I_V in iCD PC16 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc16,
          main = "TNF I_MKI vs I_V\nin iCD PC16 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_I_MKI vs I_V in iCD PC16 Genes.png')

DEG from PC1, PC14 and PC16. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.

pc1_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p1,]
pc14_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p14,]
pc16_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p16,]

pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1  & abs(pc1_dif$log2FoldChange > 0.1))
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1  & abs(pc14_dif$log2FoldChange > 0.1))
pc16_DEG <- subset(pc16_dif, pc16_dif$padj < 0.1  & abs(pc16_dif$log2FoldChange > 0.1))

pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]
pc16_genes_DEG <- pc16_genes[pc16_genes$X1 %in% pc16_DEG$X1,]


pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc16_genes_DEG <- inner_join(pc16_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))


print(paste("The number of MK2 inhibitor targets from TNF model from iCD PC1:", dim(pc1_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from iCD PC1: 1"
print(paste("The number of MK2 inhibitor targets from TNF model from iCD PC14:", dim(pc14_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from iCD PC14: 0"
print(paste("The number of MK2 inhibitor targets from TNF model from iCD PC16:", dim(pc16_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from iCD PC16: 0"

Let’s take a look at the only MK2i target from TNF model in iCD PC1 space.

pc1_DEG
## # A tibble: 1 x 7
##   X1                 baseMean log2FoldChange lfcSE  stat    pvalue   padj
##   <chr>                 <dbl>          <dbl> <dbl> <dbl>     <dbl>  <dbl>
## 1 ENSMUSG00000030064     347.          0.818 0.193  4.33 0.0000146 0.0113

I would also like to plot the Loading plot for this gene in PC1/PC14 space and PC1/PC16 space.

loading_cord <- iCD_PC[iCD_PC$`Ensembl Gene ID` %in% pc1_genes_DEG$X1,]
loading_cord <- loading_cord[, c('iCD_PC1', 'iCD_PC14', 'iCD_PC16', 'Ensembl Gene ID')]
loading_cord <- inner_join(loading_cord, ensembl_id[,c(2,4)], by = c("Ensembl Gene ID"="Ensembl Gene ID"))
loading_cord$`Ensembl Gene ID` = NULL
loading_cord[2,] <- loading_cord[1,]
loading_cord <- as.data.frame(loading_cord)
loading_cord[1,] <- c(0,0,0,NA)
loading_cord$iCD_PC1 <- as.numeric(loading_cord$iCD_PC1)
loading_cord$iCD_PC14 <- as.numeric(loading_cord$iCD_PC14)
loading_cord$iCD_PC16 <- as.numeric(loading_cord$iCD_PC16)

ggplot(loading_cord, aes(x=iCD_PC1, y=iCD_PC14)) + 
  geom_line(color = "red",
            arrow = arrow(ends = "first")) +
  xlim(-0.1, 0.1) +
  ylim(-0.1, 0.1) +
  ggtitle("Loading plot for TNF MK2i target in iCD PC1 and PC14 space") +
  geom_dl(aes(label = loading_cord$`Marker Symbol`), method = list(dl.trans(x = x - 1), dl.trans(y = y - 0.5), "last.points", cex = 0.8))
## Warning: Use of `loading_cord$`Marker Symbol`` is discouraged. Use `Marker
## Symbol` instead.

ggplot(loading_cord, aes(x=iCD_PC1, y=iCD_PC16)) + 
  geom_line(color = "red",
            arrow = arrow(ends = "first")) +
  xlim(-0.1, 0.1) +
  ylim(-0.1, 0.1) +
  ggtitle("Loading plot for TNF MK2i target in iCD PC1 and PC16 space") +
  geom_dl(aes(label = loading_cord$`Marker Symbol`), method = list(dl.trans(x = x - 1), dl.trans(y = y - 0.5), "last.points", cex = 0.8))
## Warning: Use of `loading_cord$`Marker Symbol`` is discouraged. Use `Marker
## Symbol` instead.

TNF model disregarding inflammation status

All genes in the list of 325 from PC1, PC14 and PC16

pc1_genes <- as.data.frame(TNF_combined_vst[TNF_combined_vst$X1 %in% ensembl_p1,])
pc14_genes <- as.data.frame(TNF_combined_vst[TNF_combined_vst$X1 %in% ensembl_p14,])
pc16_genes <-as.data.frame(TNF_combined_vst[TNF_combined_vst$X1 %in% ensembl_p16,])


for (i in 1:dim(pc1_genes)[1]) {
  pc1_genes[i,2:22] <- zscore(as.numeric(pc1_genes[i,2:22]))
}

for (i in 1:dim(pc14_genes)[1]) {
  pc14_genes[i,2:22] <- zscore(as.numeric(pc14_genes[i,2:22]))
}

for (i in 1:dim(pc16_genes)[1]) {
  pc16_genes[i,2:22] <- zscore(as.numeric(pc16_genes[i,2:22]))
}


heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:22])
heatmap_matrix_pc14 <- as.matrix(pc14_genes[,2:22])
heatmap_matrix_pc16 <- as.matrix(pc16_genes[,2:22])


png('TNF_MKI vs V in iCD PC1 Genes.png',
    width = 800,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
          main = "TNF MKI vs V\nin iCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_MKI vs V in iCD PC1 Genes.pdf',
    width = 8,
    height = 10)
heatmap.2(heatmap_matrix_pc1,
          main = "TNF MKI vs V\nin iCD PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_MKI vs V in iCD PC1 Genes.png')

png('TNF_MKI vs V in iCD PC14 Genes.png',
    width = 800,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc14,
          main = "TNF MKI vs V\nin iCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_MKI vs V in iCD PC14 Genes.pdf',
    width = 8,
    height = 10)
heatmap.2(heatmap_matrix_pc14,
          main = "TNF MKI vs V\nin iCD PC14 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_MKI vs V in iCD PC14 Genes.png')

png('TNF_MKI vs V in iCD PC16 Genes.png',
    width = 800,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc16,
          main = "TNF MKI vs V\nin iCD PC16 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_MKI vs V in iCD PC16 Genes.pdf',
    width = 8,
    height = 10)
heatmap.2(heatmap_matrix_pc16,
          main = "TNF MKI vs V\nin iCD PC16 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_MKI vs V in iCD PC16 Genes.png')

DEG from PC1, PC14 and PC16. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.

pc1_dif <- TNF_combined_dif[TNF_combined_dif$X1 %in% ensembl_p1,]
pc14_dif <- TNF_combined_dif[TNF_combined_dif$X1 %in% ensembl_p14,]
pc16_dif <- TNF_combined_dif[TNF_combined_dif$X1 %in% ensembl_p16,]


pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1  & abs(pc1_dif$log2FoldChange > 0.1))
pc14_DEG <- subset(pc14_dif, pc14_dif$padj < 0.1  & abs(pc14_dif$log2FoldChange > 0.1))
pc16_DEG <- subset(pc16_dif, pc16_dif$padj < 0.1  & abs(pc16_dif$log2FoldChange > 0.1))


pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc14_genes_DEG <- pc14_genes[pc14_genes$X1 %in% pc14_DEG$X1,]
pc16_genes_DEG <- pc16_genes[pc16_genes$X1 %in% pc16_DEG$X1,]


pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc14_genes_DEG <- inner_join(pc14_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc16_genes_DEG <- inner_join(pc16_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))


print(paste("The number of MK2 inhibitor targets from TNF model from iCD PC1:", dim(pc1_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from iCD PC1: 0"
print(paste("The number of MK2 inhibitor targets from TNF model from iCD PC14:", dim(pc14_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from iCD PC14: 1"
print(paste("The number of MK2 inhibitor targets from TNF model from iCD PC16:", dim(pc16_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from iCD PC16: 0"

Let’s take a look at the only MK2i target from TNF model in iCD PC1 space.

pc14_DEG
## # A tibble: 1 x 7
##   X1                 baseMean log2FoldChange lfcSE  stat    pvalue   padj
##   <chr>                 <dbl>          <dbl> <dbl> <dbl>     <dbl>  <dbl>
## 1 ENSMUSG00000089824     248.          0.892 0.201  4.27 0.0000200 0.0212

I would also like to plot the Loading plot for this gene in PC1/PC14 space and PC1/PC16 space.

loading_cord <- iCD_PC[iCD_PC$`Ensembl Gene ID` %in% pc14_genes_DEG$X1,]
loading_cord <- loading_cord[, c('iCD_PC1', 'iCD_PC14', 'iCD_PC16', 'Ensembl Gene ID')]
loading_cord <- inner_join(loading_cord, ensembl_id[,c(2,4)], by = c("Ensembl Gene ID"="Ensembl Gene ID"))
loading_cord$`Ensembl Gene ID` = NULL
loading_cord[2,] <- loading_cord[1,]
loading_cord <- as.data.frame(loading_cord)
loading_cord[1,] <- c(0,0,0,NA)
loading_cord$iCD_PC1 <- as.numeric(loading_cord$iCD_PC1)
loading_cord$iCD_PC14 <- as.numeric(loading_cord$iCD_PC14)
loading_cord$iCD_PC16 <- as.numeric(loading_cord$iCD_PC16)

ggplot(loading_cord, aes(x=iCD_PC1, y=iCD_PC14)) + 
  geom_line(color = "red",
            arrow = arrow(ends = "first")) +
  xlim(-0.1, 0.1) +
  ylim(-0.1, 0.1) +
  ggtitle("Loading plot for TNF MK2i target in iCD PC1 and PC14 space") +
  geom_dl(aes(label = loading_cord$`Marker Symbol`), method = list(dl.trans(x = x - 1), dl.trans(y = y - 0.5), "last.points", cex = 0.8))
## Warning: Use of `loading_cord$`Marker Symbol`` is discouraged. Use `Marker
## Symbol` instead.

ggplot(loading_cord, aes(x=iCD_PC1, y=iCD_PC16)) + 
  geom_line(color = "red",
            arrow = arrow(ends = "first")) +
  xlim(-0.1, 0.1) +
  ylim(-0.1, 0.1) +
  ggtitle("Loading plot for TNF MK2i target in iCD PC1 and PC16 space") +
  geom_dl(aes(label = loading_cord$`Marker Symbol`), method = list(dl.trans(x = x - 1), dl.trans(y = y - 0.5), "last.points", cex = 0.8))
## Warning: Use of `loading_cord$`Marker Symbol`` is discouraged. Use `Marker
## Symbol` instead.

Heatmap plotting based on UC PC using MK2i dataset

TCT and TNF MK2i treated dataset with and without inflammation are used here.

# ID conversion
# load the PC dataset
UC_PC <- read_csv("UC PC.csv")
## Parsed with column specification:
## cols(
##   Row = col_character(),
##   Var1_1 = col_double(),
##   Var1_2 = col_double(),
##   Var1_3 = col_double(),
##   Var1_4 = col_double(),
##   Var1_5 = col_double(),
##   Var1_6 = col_double(),
##   Var1_7 = col_double(),
##   Var1_8 = col_double(),
##   Var1_9 = col_double(),
##   Var1_10 = col_double()
## )
human_id <- as.data.frame(UC_PC$Row)
human_id <- as_tibble(human_id)
colnames(human_id) <- "human_gene_symbol"

id_conversion <- inner_join(human_id, homolog_id, by = c("human_gene_symbol" = "Human Marker Symbol"))
## Warning: Column `human_gene_symbol`/`Human Marker Symbol` joining factor and
## character vector, coercing into character vector
non_duplicates_idx <- which(duplicated(id_conversion$human_gene_symbol) == FALSE)
id_conversion <- id_conversion[non_duplicates_idx, ]

id_conversion <- inner_join(id_conversion, ensembl_id, by = c("Mouse Marker Symbol" = "Marker Symbol"))
non_duplicates_idx <- which(duplicated(id_conversion$`Mouse Marker Symbol`) == FALSE)
id_conversion <- id_conversion[non_duplicates_idx, ]

id_conversion <- id_conversion[,c(1,7)]

UC_PC <- inner_join(UC_PC, id_conversion, by = c("Row" = "human_gene_symbol"))

For UC PCs, PC1 and PC2 are selected because they seperates MK2i treatment statuses in TCT mice.

The goal of this is to examine how many of the genes that contribute to the seperation of MK2i vs Vehicle treated TCT mice based on human PCs are targeted by MK2 inhibitor.

For each PC, the top 49 contributing genes are selected for plotting the heatmap (10% of the list). Top contributing genes are genes with the largest absolute loading coefficients.

n = 10
threshold_p1 <- quantile(abs(UC_PC$Var1_1), prob=1-n/100)
threshold_p2 <- quantile(abs(UC_PC$Var1_2), prob=1-n/100)

ensembl_p1 <- UC_PC$`Ensembl Gene ID`[abs(UC_PC$Var1_1) > threshold_p1]
ensembl_p2 <- UC_PC$`Ensembl Gene ID`[abs(UC_PC$Var1_2) > threshold_p2]

TCT model with inflammation treated with MK2i

All genes in the list of 49 from PC1 and PC2

pc1_genes <- as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p1,])
pc2_genes <- as.data.frame(TCT_imki_iv_vst[TCT_imki_iv_vst$X1 %in% ensembl_p2,])


for (i in 1:dim(pc1_genes)[1]) {
  pc1_genes[i,2:13] <- zscore(as.numeric(pc1_genes[i,2:13]))
}

for (i in 1:dim(pc2_genes)[1]) {
  pc2_genes[i,2:13] <- zscore(as.numeric(pc2_genes[i,2:13]))
}

heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:13])
heatmap_matrix_pc2 <- as.matrix(pc2_genes[,2:13])


png('TCT_I_MKI vs I_V in UC PC1 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
          main = "TCT I_MKI vs I_V\nin UC PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_I_MKI vs I_V in UC PC1 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc1,
          main = "TCT I_MKI vs I_V\nin UC PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TCT_I_MKI vs I_V in UC PC1 Genes.png')

png('TCT_I_MKI vs I_V in UC PC2 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc2,
          main = "TCT I_MKI vs I_V\nin UC PC2 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_I_MKI vs I_V in UC PC2 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc2,
          main = "TCT I_MKI vs I_V\nin UC PC2 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TCT_I_MKI vs I_V in UC PC2 Genes.png')

DEG from PC1 and PC2. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.

pc1_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p1,]
pc2_dif <- TCT_imki_iv_dif[TCT_imki_iv_dif$X1 %in% ensembl_p2,]
pc1_dif
## # A tibble: 49 x 7
##    X1                 baseMean log2FoldChange lfcSE   stat pvalue  padj
##    <chr>                 <dbl>          <dbl> <dbl>  <dbl>  <dbl> <dbl>
##  1 ENSMUSG00000004994     84.8         -0.508 0.295 -1.57   0.115  1.00
##  2 ENSMUSG00000005299   1920.          -0.119 0.184 -0.587  0.557  1.00
##  3 ENSMUSG00000006205    486.           0.318 0.257  1.09   0.274  1.00
##  4 ENSMUSG00000020212     64.7         -0.344 0.268 -1.22   0.223  1.00
##  5 ENSMUSG00000020541    839.          -0.129 0.188 -0.613  0.540  1.00
##  6 ENSMUSG00000020719   8112.          -0.220 0.247 -0.797  0.425  1.00
##  7 ENSMUSG00000020810    291.           0.194 0.195  0.883  0.377  1.00
##  8 ENSMUSG00000020865   1996.           0.482 0.263  1.52   0.128  1.00
##  9 ENSMUSG00000021356     44.2          0.343 0.290  0.879  0.379  1.00
## 10 ENSMUSG00000021428    224.          -0.313 0.233 -1.19   0.234  1.00
## # … with 39 more rows
pc2_dif
## # A tibble: 49 x 7
##    X1                 baseMean log2FoldChange lfcSE    stat pvalue  padj
##    <chr>                 <dbl>          <dbl> <dbl>   <dbl>  <dbl> <dbl>
##  1 ENSMUSG00000000441    1301.        -0.364  0.147 -2.44   0.0148  1.00
##  2 ENSMUSG00000001870     236.         0.0587 0.232  0.0630 0.950   1.00
##  3 ENSMUSG00000008575     628.         0.262  0.236  0.918  0.359   1.00
##  4 ENSMUSG00000009555    1008.        -0.139  0.169 -0.777  0.437   1.00
##  5 ENSMUSG00000010392     766.        -0.308  0.244 -1.09   0.275   1.00
##  6 ENSMUSG00000012350    3401.        -0.172  0.259 -0.528  0.597   1.00
##  7 ENSMUSG00000017418     579.        -0.481  0.243 -1.87   0.0617  1.00
##  8 ENSMUSG00000017478     871.        -0.351  0.278 -1.09   0.275   1.00
##  9 ENSMUSG00000019960    1159.        -0.336  0.240 -1.30   0.192   1.00
## 10 ENSMUSG00000020149    5883.        -0.163  0.131 -1.23   0.220   1.00
## # … with 39 more rows
pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1  & abs(pc1_dif$log2FoldChange > 0.1))
pc2_DEG <- subset(pc2_dif, pc2_dif$padj < 0.1  & abs(pc2_dif$log2FoldChange > 0.1))

pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc2_genes_DEG <- pc2_genes[pc2_genes$X1 %in% pc2_DEG$X1,]

pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc2_genes_DEG <- inner_join(pc2_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))

print(paste("The number of MK2 inhibitor targets from TCT model from UC PC1:", dim(pc1_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TCT model from UC PC1: 0"
print(paste("The number of MK2 inhibitor targets from TCT model from UC PC2:", dim(pc2_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TCT model from UC PC2: 0"

TCT model with no inflammation treated with MK2i

All genes in the list of 49 from PC1 and PC2

pc1_genes <- as.data.frame(TCT_nimki_niv_vst[TCT_nimki_niv_vst$X1 %in% ensembl_p1,])
pc2_genes <- as.data.frame(TCT_nimki_niv_vst[TCT_nimki_niv_vst$X1 %in% ensembl_p2,])


for (i in 1:dim(pc1_genes)[1]) {
  pc1_genes[i,2:12] <- zscore(as.numeric(pc1_genes[i,2:12]))
}

for (i in 1:dim(pc2_genes)[1]) {
  pc2_genes[i,2:12] <- zscore(as.numeric(pc2_genes[i,2:12]))
}

heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:12])
heatmap_matrix_pc2 <- as.matrix(pc2_genes[,2:12])


png('TCT_NI_MKI vs NI_V in UC PC1 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
          main = "TCT NI_MKI vs NI_V\nin UC PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_NI_MKI vs NI_V in UC PC1 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc1,
          main = "TCT NI_MKI vs NI_V\nin UC PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TCT_NI_MKI vs NI_V in UC PC1 Genes.png')

png('TCT_NI_MKI vs NI_V in UC PC2 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc2,
          main = "TCT NI_MKI vs NI_V\nin UC PC2 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TCT_NI_MKI vs NI_V in UC PC2 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc2,
          main = "TCT NI_MKI vs NI_V\nin UC PC2 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TCT_NI_MKI vs NI_V in UC PC2 Genes.png')

DEG from PC1 and PC2. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.

pc1_dif <- TCT_nimki_niv_dif[TCT_nimki_niv_dif$X1 %in% ensembl_p1,]
pc2_dif <- TCT_nimki_niv_dif[TCT_nimki_niv_dif$X1 %in% ensembl_p2,]
pc1_dif
## # A tibble: 49 x 7
##    X1                 baseMean log2FoldChange lfcSE    stat pvalue  padj
##    <chr>                 <dbl>          <dbl> <dbl>   <dbl>  <dbl> <dbl>
##  1 ENSMUSG00000004994    102.         -0.296  0.191 -1.65   0.0981  1.00
##  2 ENSMUSG00000005299   2443.         -0.107  0.126 -0.823  0.410   1.00
##  3 ENSMUSG00000006205    507.          0.0214 0.188  0.0193 0.985   1.00
##  4 ENSMUSG00000020212     54.0         0.0328 0.194  0.118  0.906   1.00
##  5 ENSMUSG00000020541   1084.          0.131  0.101  1.32   0.186   1.00
##  6 ENSMUSG00000020719   9771.         -0.0156 0.146 -0.0439 0.965   1.00
##  7 ENSMUSG00000020810    255.         -0.144  0.166 -0.877  0.381   1.00
##  8 ENSMUSG00000020865   3098.         -0.126  0.144 -0.892  0.372   1.00
##  9 ENSMUSG00000021356     12.2        -0.0571 0.172 -0.559  0.576   1.00
## 10 ENSMUSG00000021428    200.          0.0337 0.156  0.243  0.808   1.00
## # … with 39 more rows
pc2_dif
## # A tibble: 49 x 7
##    X1                 baseMean log2FoldChange  lfcSE    stat pvalue  padj
##    <chr>                 <dbl>          <dbl>  <dbl>   <dbl>  <dbl> <dbl>
##  1 ENSMUSG00000000441    1486.        0.0277  0.142   0.159   0.873  1.00
##  2 ENSMUSG00000001870     319.        0.0423  0.161   0.183   0.855  1.00
##  3 ENSMUSG00000008575    1119.       -0.0938  0.115  -0.768   0.443  1.00
##  4 ENSMUSG00000009555     999.        0.0492  0.142   0.334   0.739  1.00
##  5 ENSMUSG00000010392    1012.       -0.0414  0.191  -0.0328  0.974  1.00
##  6 ENSMUSG00000012350    6681.       -0.00547 0.165   0.0936  0.925  1.00
##  7 ENSMUSG00000017418     706.       -0.135   0.175  -0.613   0.540  1.00
##  8 ENSMUSG00000017478    1033.        0.133   0.176   0.851   0.395  1.00
##  9 ENSMUSG00000019960     776.       -0.00514 0.193  -0.0380  0.970  1.00
## 10 ENSMUSG00000020149    6450.       -0.120   0.0941 -1.26    0.207  1.00
## # … with 39 more rows
pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1  & abs(pc1_dif$log2FoldChange > 0.1))
pc2_DEG <- subset(pc2_dif, pc2_dif$padj < 0.1  & abs(pc2_dif$log2FoldChange > 0.1))

pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc2_genes_DEG <- pc2_genes[pc2_genes$X1 %in% pc2_DEG$X1,]

pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc2_genes_DEG <- inner_join(pc2_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))

print(paste("The number of MK2 inhibitor targets from TCT model from UC PC1:", dim(pc1_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TCT model from UC PC1: 0"
print(paste("The number of MK2 inhibitor targets from TCT model from UC PC2:", dim(pc2_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TCT model from UC PC2: 0"

TNF model with inflammation treated with MK2i

All genes in the list of 49 from PC1 and PC2

pc1_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p1,])
pc2_genes <- as.data.frame(TNF_imki_iv_vst[TNF_imki_iv_vst$X1 %in% ensembl_p2,])


for (i in 1:dim(pc1_genes)[1]) {
  pc1_genes[i,2:12] <- zscore(as.numeric(pc1_genes[i,2:12]))
}

for (i in 1:dim(pc2_genes)[1]) {
  pc2_genes[i,2:12] <- zscore(as.numeric(pc2_genes[i,2:12]))
}

heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:12])
heatmap_matrix_pc2 <- as.matrix(pc2_genes[,2:12])


png('TNF_I_MKI vs I_V in UC PC1 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
          main = "TNF I_MKI vs I_V\nin UC PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_I_MKI vs I_V in UC PC1 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc1,
          main = "TNF I_MKI vs I_V\nin UC PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_I_MKI vs I_V in UC PC1 Genes.png')

png('TNF_I_MKI vs I_V in UC PC2 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc2,
          main = "TNF I_MKI vs I_V\nin UC PC2 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_I_MKI vs I_V in UC PC2 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc2,
          main = "TNF I_MKI vs I_V\nin UC PC2 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_I_MKI vs I_V in UC PC2 Genes.png')

DEG from PC1 and PC2. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.

pc1_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p1,]
pc2_dif <- TNF_imki_iv_dif[TNF_imki_iv_dif$X1 %in% ensembl_p2,]
pc1_dif
## # A tibble: 49 x 7
##    X1                 baseMean log2FoldChange lfcSE   stat pvalue  padj
##    <chr>                 <dbl>          <dbl> <dbl>  <dbl>  <dbl> <dbl>
##  1 ENSMUSG00000004994     48.9         0.150  0.246  0.603 0.547   1.00
##  2 ENSMUSG00000005299   1109.         -0.126  0.108 -1.14  0.255   1.00
##  3 ENSMUSG00000006205    470.          0.0540 0.240 -0.148 0.883   1.00
##  4 ENSMUSG00000020212     31.5        -0.0306 0.239 -0.217 0.828   1.00
##  5 ENSMUSG00000020541    503.         -0.0682 0.157 -0.496 0.620   1.00
##  6 ENSMUSG00000020719   5321.         -0.216  0.119 -1.79  0.0738  1.00
##  7 ENSMUSG00000020810    148.          0.366  0.212  1.62  0.105   1.00
##  8 ENSMUSG00000020865    903.         -0.0821 0.168 -0.339 0.735   1.00
##  9 ENSMUSG00000021356    242.          0.162  0.203  0.558 0.577   1.00
## 10 ENSMUSG00000021428    119.         -0.111  0.187 -0.455 0.649   1.00
## # … with 39 more rows
pc2_dif
## # A tibble: 49 x 7
##    X1                 baseMean log2FoldChange lfcSE    stat pvalue  padj
##    <chr>                 <dbl>          <dbl> <dbl>   <dbl>  <dbl> <dbl>
##  1 ENSMUSG00000000441     751.        -0.164  0.142 -1.08   0.278   1.00
##  2 ENSMUSG00000001870     232.         0.375  0.235  1.22   0.224   1.00
##  3 ENSMUSG00000008575     442.         0.0164 0.202  0.0658 0.948   1.00
##  4 ENSMUSG00000009555     575.         0.264  0.184  1.47   0.140   1.00
##  5 ENSMUSG00000010392     499.         0.0190 0.151 -0.122  0.903   1.00
##  6 ENSMUSG00000012350     967.         0.382  0.209  1.66   0.0960  1.00
##  7 ENSMUSG00000017418     311.         0.0123 0.207  0.0503 0.960   1.00
##  8 ENSMUSG00000017478     538.         0.210  0.242  0.949  0.342   1.00
##  9 ENSMUSG00000019960     584.        -0.287  0.193 -1.32   0.187   1.00
## 10 ENSMUSG00000020149    3198.        -0.127  0.153 -0.873  0.383   1.00
## # … with 39 more rows
pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1  & abs(pc1_dif$log2FoldChange > 0.1))
pc2_DEG <- subset(pc2_dif, pc2_dif$padj < 0.1  & abs(pc2_dif$log2FoldChange > 0.1))

pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc2_genes_DEG <- pc2_genes[pc2_genes$X1 %in% pc2_DEG$X1,]

pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc2_genes_DEG <- inner_join(pc2_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))

print(paste("The number of MK2 inhibitor targets from TNF model from UC PC1:", dim(pc1_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from UC PC1: 0"
print(paste("The number of MK2 inhibitor targets from TNF model from UC PC2:", dim(pc2_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from UC PC2: 0"

TNF model with no inflammation treated with MK2i

All genes in the list of 49 from PC1 and PC2

pc1_genes <- as.data.frame(TNF_nimki_niv_vst[TNF_nimki_niv_vst$X1 %in% ensembl_p1,])
pc2_genes <- as.data.frame(TNF_nimki_niv_vst[TNF_nimki_niv_vst$X1 %in% ensembl_p2,])


for (i in 1:dim(pc1_genes)[1]) {
  pc1_genes[i,2:11] <- zscore(as.numeric(pc1_genes[i,2:11]))
}

for (i in 1:dim(pc2_genes)[1]) {
  pc2_genes[i,2:11] <- zscore(as.numeric(pc2_genes[i,2:11]))
}

heatmap_matrix_pc1 <- as.matrix(pc1_genes[,2:11])
heatmap_matrix_pc2 <- as.matrix(pc2_genes[,2:11])


png('TNF_NI_MKI vs NI_V in UC PC1 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc1,
          main = "TNF NI_MKI vs NI_V\nin UC PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_NI_MKI vs NI_V in UC PC1 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc1,
          main = "TNF NI_MKI vs NI_V\nin UC PC1 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_NI_MKI vs NI_V in UC PC1 Genes.png')

png('TNF_NI_MKI vs NI_V in UC PC2 Genes.png',
    width = 600,
    height = 1200,
    res = 200,
    pointsize = 8)
heatmap.2(heatmap_matrix_pc2,
          main = "TNF NI_MKI vs NI_V\nin UC PC2 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
pdf('TNF_NI_MKI vs NI_V in UC PC2 Genes.pdf',
    width = 6,
    height = 10)
heatmap.2(heatmap_matrix_pc2,
          main = "TNF NI_MKI vs NI_V\nin UC PC2 Genes",
          density.info = "none",
          key = TRUE,
          lhei = c(1,7),
          col=my_palette,
          cexCol = 1,
          margins = c(8,2),
          trace = "none",
          dendrogram = "row",
          labRow = FALSE,
          keysize = 2,
          ylab = "Genes",
          Colv = "NA")
dev.off()
## quartz_off_screen 
##                 2
include_graphics('TNF_NI_MKI vs NI_V in UC PC2 Genes.png')

DEG from PC1 and PC2. The standard of DEG is relaxed to p < 0.1 here to keep consistent with previous analysis and LFC > 0.1.

pc1_dif <- TNF_nimki_niv_dif[TNF_nimki_niv_dif$X1 %in% ensembl_p1,]
pc2_dif <- TNF_nimki_niv_dif[TNF_nimki_niv_dif$X1 %in% ensembl_p2,]
pc1_dif
## # A tibble: 49 x 7
##    X1                 baseMean log2FoldChange lfcSE    stat pvalue  padj
##    <chr>                 <dbl>          <dbl> <dbl>   <dbl>  <dbl> <dbl>
##  1 ENSMUSG00000004994     73.1        0.0694  0.198  0.216   0.829  1.00
##  2 ENSMUSG00000005299   1924.         0.00207 0.146  0.0132  0.989  1.00
##  3 ENSMUSG00000006205    202.        -0.294   0.198 -1.43    0.153  1.00
##  4 ENSMUSG00000020212     33.8        0.134   0.192  0.633   0.527  1.00
##  5 ENSMUSG00000020541    745.         0.0694  0.148  0.449   0.654  1.00
##  6 ENSMUSG00000020719   6565.        -0.00388 0.126 -0.0323  0.974  1.00
##  7 ENSMUSG00000020810    125.        -0.124   0.193 -0.534   0.593  1.00
##  8 ENSMUSG00000020865    681.        -0.0432  0.186 -0.230   0.818  1.00
##  9 ENSMUSG00000021356     99.6       -0.133   0.193 -0.744   0.457  1.00
## 10 ENSMUSG00000021428    142.        -0.0579  0.156 -0.372   0.710  1.00
## # … with 39 more rows
pc2_dif
## # A tibble: 49 x 7
##    X1                 baseMean log2FoldChange lfcSE    stat pvalue  padj
##    <chr>                 <dbl>          <dbl> <dbl>   <dbl>  <dbl> <dbl>
##  1 ENSMUSG00000000441    1066.       -0.0407  0.155 -0.269   0.788  1.00
##  2 ENSMUSG00000001870     179.        0.0645  0.198  0.552   0.581  1.00
##  3 ENSMUSG00000008575     523.       -0.158   0.196 -0.725   0.468  1.00
##  4 ENSMUSG00000009555     612.        0.216   0.165  1.29    0.198  1.00
##  5 ENSMUSG00000010392     679.        0.00609 0.144  0.0390  0.969  1.00
##  6 ENSMUSG00000012350    1293.        0.106   0.195  0.579   0.563  1.00
##  7 ENSMUSG00000017418     377.       -0.149   0.183 -0.827   0.408  1.00
##  8 ENSMUSG00000017478     698.        0.0420  0.196  0.227   0.821  1.00
##  9 ENSMUSG00000019960     736.        0.0174  0.185  0.0831  0.934  1.00
## 10 ENSMUSG00000020149    4293.       -0.0580  0.112 -0.531   0.596  1.00
## # … with 39 more rows
pc1_DEG <- subset(pc1_dif, pc1_dif$padj < 0.1  & abs(pc1_dif$log2FoldChange > 0.1))
pc2_DEG <- subset(pc2_dif, pc2_dif$padj < 0.1  & abs(pc2_dif$log2FoldChange > 0.1))

pc1_genes_DEG <- pc1_genes[pc1_genes$X1 %in% pc1_DEG$X1,]
pc2_genes_DEG <- pc2_genes[pc2_genes$X1 %in% pc2_DEG$X1,]

pc1_genes_DEG <- inner_join(pc1_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))
pc2_genes_DEG <- inner_join(pc2_genes_DEG, ensembl_id[,c(2,4)], by = c("X1" = "Ensembl Gene ID"))

print(paste("The number of MK2 inhibitor targets from TNF model from UC PC1:", dim(pc1_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from UC PC1: 0"
print(paste("The number of MK2 inhibitor targets from TNF model from UC PC2:", dim(pc2_genes_DEG)[1]))
## [1] "The number of MK2 inhibitor targets from TNF model from UC PC2: 0"